library(KFAS)
library(xts)
library(forecast)
library(Metrics)
library(reticulate)
use_python("/usr/local/python")
Caricamento della serie storica e conversione dei dati in oggetti xts
dt <- read.csv("csv_weights/time_series_dataset.csv", sep = ";")
y <- dt$value
y <- xts(y, as.Date(as.character(dt$Data), format = "%Y-%m-%d"))
plot(y)
Confronto media e stdev per valutare la trasformazione logaritmica per sistemare la stazionarietà in varianza.
means <- as.numeric(apply.monthly(y, mean))
stdev <- as.numeric(apply.monthly(y, sd))
plot(means, stdev)
abline(lm(stdev ~ means), col = "red")
Trasformazione di BoxCox per risolvere la non stazionarietà in varianza
y_box <- BoxCox(y, lambda = "auto")
lambda <- BoxCox.lambda(y)
y_train <- y_box["2010-01-01/2017-12-31"]
y_train_num <- as.numeric(y_train)
y_val <- y_box["2018-01-01/2018-12-31"]
y_val_num <- as.numeric(y_val)
#write.csv(y_box, "time_series_trasformed.csv")
ggtsdisplay(y_train, lag.max = 100)
ggtsdisplay(y_train, lag.max = 1000)
I grafici vengono plottati con due diversi lag.max in modo da cogliere sia gli andamenti a breve termine che a lungo termine. Da acf può notare una discesa lineare e la presenza di stagionalità settimanale e annua Da pacf si può notare una discesa geometrica ogni 7 giorni Guardando i grafici si può cogliere un trend che rende la serie non stazionaria in media. Il trend sembra essere inizialmente descrescente e poi crescente. Eseguiamo le due integrazioni:
mod1 <- Arima(y_train, order = c(0,1,0), seasonal = list(order = c(0,1,0), period= 7))
ggtsdisplay(mod1$residuals, lag.max = 50)
ggtsdisplay(mod1$residuals, lag.max = 1000)
Guardando la Pacf sembra esserci sia una componente AR(p)MA(p) che una componente SMA(7)
mod2 <- Arima(y_train, order = c(1,1,1), seasonal = list(order = c(0,1,1), period= 7))
ggtsdisplay(mod2$residuals, lag.max = 100)
#ggtsdisplay(mod2$residuals, lag.max = 1000)
I grafici suggeriscono l’inserimento di un SAR(7)
mod3 <- Arima(y_train, order = c(1,1,1), seasonal = list(order = c(1,1,1), period= 7))
checkresiduals(mod3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,1,1)[7]
## Q* = 21.54, df = 6, p-value = 0.001467
##
## Model df: 4. Total lags used: 10
ggtsdisplay(mod3$residuals, lag.max = 100)
ggtsdisplay(mod3$residuals, lag.max = 1000)
I lag sembrano essere sotto le bande al 95%. Rimane da gestirle la stagionalità annua e si può effettuare una sorta di grid search per aumentare la bontà del modello che tenga in considerazione il numero di termini di fourier per cogliere l’informazione intra-annua e l’ordine di AR. La scelta del modell migliore viene fatta guardando l’AIC.
i <- 1
harm <- c(1,5,10)
best_arima <- data.frame(p = rep(NA,18), k = rep(NA,18), AIC = rep(NA,18))
for(p in seq(1:6)) {
for (k in harm) {
fs_train <- fourier(ts(y_train, frequency = 365), K = k)
fs_test <- fourier(ts(y_train, frequency = 365), K = k, h = 365)
an.error.occured <- FALSE
tryCatch( { model <- Arima(y_train, c(p,1,1), list(order=c(1,1,1), period=7),
xreg=fs_train, include.constant = TRUE); print(model$aic) },
error = function(e) {an.error.occured <<- TRUE})
if(an.error.occured == FALSE) {
best_arima[i,] <- c(p, k, model$aic)
}
i <- i + 1
}
}
best_arima
write.csv(best_arima,"csv_weights/best_arima.csv", row.names = FALSE)
best_arima <- read.csv("csv_weights/best_arima.csv")
best_arima_config <- best_arima[which.min(best_arima$AIC),]
best_arima_config
## p k AIC
## 18 6 10 22012.17
fs_train <- fourier(ts(y_train, frequency = 365), K = best_arima_config$k)
fs_test <- fourier(ts(y_train, frequency = 365), K = best_arima_config$k, h = 365)
mod4 <- Arima(y_train, c(best_arima_config$p,1,1), list(order=c(1,1,1), period=7),
xreg=fs_train, include.constant = TRUE)
pred4<- forecast(mod4, h=365, xreg=fs_test)
autoplot(pred4)
print(paste0("MAPE on training set ", mean(abs(mod4$fitted - y_train_num)/y_train_num) * 100))
## [1] "MAPE on training set 8.99792859533606"
print(paste0("MAPE on validation set ", mean(abs(pred4$mean - y_val_num)/y_val_num) * 100))
## [1] "MAPE on validation set 10.9370875875956"
plot(y_val_num, type = "l",)
lines(as.numeric(pred4$mean), col = "blue")
Inserimento delle festività tramite dummy:
library(dsa)
dm <- data.frame(Date = seq(from = as.Date("2010-01-01"), to = as.Date("2019-12-31"), by = 1))
holidays <- c("EasterMon","EasterTue","EasterSat","Nov1","Dec24","Dec25","Dec26","Dec31","Jan1","Jan6","Dec8","Aug15","Apr25","Mag1","Giu2")
#timeDate::listHolidays()
EasterMon <- as.Date(c(timeDate::EasterMonday(2010:2019)))
EasterTue <- EasterMon + 1
EasterSat <- EasterMon - 2
Nov1 <- as.Date(c(timeDate::AllSaints(2010:2019)))
Dec24 <- as.Date(c(timeDate::ChristmasEve(2010:2019)))
Dec25 <- Dec24 + 1
Dec26 <- Dec25 + 1
Dec31 <- Dec26 + 5
Jan1 <- Dec26 + 6 - 365
Jan6 <- Jan1 + 5
Dec8 <- as.Date(c(timeDate::ITImmaculateConception(2010:2019)))
Aug15 <- as.Date(c(timeDate::ITAssumptionOfVirginMary(2010:2019)))
Apr25 <- as.Date(c(timeDate::ITLiberationDay(2010:2019)))
Mag1 <- Apr25 + 6
Giu2 <- Mag1 + 32
for( i in holidays) {
holy <- rep(0,length(dm$Date))
dm <- cbind(dm, holy)
}
colnames(dm) <- c("Date",holidays)
for( j in seq(10)) {
dm$EasterTue[which(dm$Date == EasterTue[j])] <- 1
dm$EasterMon[which(dm$Date == EasterMon[j])] <- 1
dm$EasterSat[which(dm$Date == EasterSat[j])] <- 1
dm$Nov1[which(dm$Date == Nov1[j])] <- 1
dm$Dec24[which(dm$Date == Dec24[j])] <- 1
dm$Dec25[which(dm$Date == Dec25[j])] <- 1
dm$Dec26[which(dm$Date == Dec26[j])] <- 1
dm$Dec31[which(dm$Date == Dec31[j])] <- 1
dm$Jan1[which(dm$Date == Jan1[j])] <- 1
dm$Jan6[which(dm$Date == Jan6[j])] <- 1
dm$Dec8[which(dm$Date == Dec8[j])] <- 1
dm$Aug15[which(dm$Date == Aug15[j])] <- 1
dm$Apr25[which(dm$Date == Apr25[j])] <- 1
dm$Mag1[which(dm$Date == Mag1[j])] <- 1
dm$Giu2[which(dm$Date == Giu2[j])] <- 1
}
hd <- rowSums(dm[holidays])
dm <- cbind(dm, hd)
dm$weekday <- weekdays(as.Date(dm$Date))
dm[which(dm$weekday == "Domenica"),][holidays] <- c(rep(0,15))
write.csv(dm,"csv_weights/dummy_project.csv", row.names = FALSE)
holidays <- c("EasterMon","EasterTue","EasterSat","Nov1","Dec24","Dec25","Dec26","Dec31","Jan1","Jan6","Dec8","Aug15","Apr25","Mag1","Giu2")
dm <-read.csv("csv_weights/dummy_project.csv")
fs_train <- fourier(ts(y_train, frequency = 365), K = best_arima_config$k)
fs_test <- fourier(ts(y_train, frequency = 365), K = best_arima_config$k, h = 365)
dummy_holidays_train <- dm[holidays][1:length(y_train),]
dummy_holidays_val <- dm[holidays][(length(y_train)+1):(length(y_train)+365),]
mod4_reg <- Arima(y_train, c(best_arima_config$p,1,1), list(order=c(1,1,1), period=7),
xreg=as.matrix(cbind(fs_train,dummy_holidays_train)), include.constant = TRUE)
#checkresiduals(mod4_reg)
pred4_reg <- forecast(mod4_reg, h=365, xreg=as.matrix(cbind(fs_test,dummy_holidays_val)))
#autoplot(pred4_reg)
print(paste0("MAPE on training set ", mean(abs(mod4_reg$fitted - y_train_num)/y_train_num) * 100))
## [1] "MAPE on training set 8.45300064992593"
print(paste0("MAPE on validation set ", mean(abs(pred4_reg$mean - y_val_num)/y_val_num) * 100))
## [1] "MAPE on validation set 10.918625828688"
# plot validation
plot(y_val_num, type = "l")
lines(as.numeric(pred4_reg$mean), col = "blue")
Proviamo a far variare k tenendo lo stesso modello e le dummy delle festivitÃ
i <- 1
harm <- c(0,1,2,4,6,10)
best_arima_reg <- data.frame(k = rep(NA,6), AIC = rep(NA,6), Mape_val = rep(NA,6))
for (k in harm) {
an.error.occured <- FALSE
if(k != 0) {
fs_train <- fourier(ts(y_train, frequency = 365), K = k)
fs_test <- fourier(ts(y_train, frequency = 365), K = k, h = 365)
tryCatch( { model <- Arima(y_train, c(best_arima_config$p,1,1), list(order=c(1,1,1), period=7),
xreg=as.matrix(cbind(fs_train,dummy_holidays_train)), include.constant = TRUE); print(model$aic) },
error = function(e) {an.error.occured <<- TRUE})
prediction <- forecast(model, h=365, xreg=as.matrix(cbind(fs_test,dummy_holidays_val)))
} else {
model <- Arima(y_train, c(best_arima_config$p,1,1), list(order=c(1,1,1), period=7),
xreg=as.matrix(dummy_holidays_train), include.constant = TRUE)
prediction <- forecast(model, h=365, xreg=as.matrix(dummy_holidays_val))
}
if(an.error.occured == FALSE) {
mape_val <- mean(abs(prediction$mean - y_val_num)/y_val_num) * 100
best_arima_reg[i,] <- c(k, model$aic, mape_val)
}
i <- i + 1
}
best_arima_reg
write.csv(best_arima_reg,"csv_weights/best_arima_reg.csv", row.names = FALSE)
best_arima_reg <- read.csv("csv_weights/best_arima_reg.csv")
best_arima_reg_config <- best_arima_reg[which.min(best_arima_reg$Mape_val),]
best_arima_reg_config
## k AIC Mape_val
## 2 1 21613.34 10.08227
Miglior modello arima con regressori dummy per la festività e K = 1 termini di Fourier
fs_train <- fourier(ts(y_train, frequency = 365), K = best_arima_reg_config$k)
fs_test <- fourier(ts(y_train, frequency = 365), K = best_arima_reg_config$k, h = 365)
dummy_holidays_train <- dm[holidays][1:length(y_train),]
dummy_holidays_val <- dm[holidays][(length(y_train)+1):(length(y_train)+365),]
mod5_reg <- Arima(y_train, c(best_arima_config$p,1,1), list(order=c(1,1,1), period=7),
xreg=as.matrix(cbind(fs_train,dummy_holidays_train)), include.constant = TRUE)
#mod5_reg
#checkresiduals(mod5_reg)
pred5_reg <- forecast(mod5_reg, h=365, xreg=as.matrix(cbind(fs_test,dummy_holidays_val)))
autoplot(pred5_reg)
print(paste0("MAPE on training set ", mean(abs(mod5_reg$fitted - y_train_num)/y_train_num) * 100))
## [1] "MAPE on training set 8.53247472403642"
print(paste0("MAPE on validation set ", mean(abs(pred5_reg$mean - y_val_num)/y_val_num) * 100))
## [1] "MAPE on validation set 10.0822667905118"
# plot totale
plot(as.numeric(y_box), type ="l")
lines(pred5_reg$fitted, col = "red")
lines(pred5_reg$mean, col = "blue")
# plot validation
plot(y_val_num, type = "l")
lines(as.numeric(pred4$mean), col = "blue")
lines(as.numeric(pred5_reg$mean), col = "red")
legend(280, 160, legend=c("whitout dummy", "with dummy"), col=c("blue", "red"), lty=1:1, cex=0.8)
checkresiduals(mod4_reg)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(6,1,1)(1,1,1)[7] errors
## Q* = 78.634, df = 3, p-value < 2.2e-16
##
## Model df: 44. Total lags used: 47
checkresiduals(mod5_reg)
##
## Ljung-Box test
##
## data: Residuals from Regression with ARIMA(6,1,1)(1,1,1)[7] errors
## Q* = 54.442, df = 3, p-value = 9.033e-12
##
## Model df: 26. Total lags used: 29
Date le caratteristiche della serie storica sono stati provati due modelli UCM: - RW con stagionalità settimanale (dummy) e annuale (trigonometrica) - IRW con stagionalità settimale (dummy) e annuale (trigonometrica)
Su entrambi i modelli è stato effettuato un grid search per scegliere il numero di armoniche migliore
y_ucm <- y_train_num
y_ucm[(length(y_train)+1):(nrow(dt))] <- NA
#plot(y_ucm, type ="l")
best_ucm <- data.frame(Name = rep(NA,8), Mape_train = rep(NA,8), Mape_val = rep(NA,8))
harm <- c(2,6,10,12,14)
vary <- var(y_ucm, na.rm = TRUE)
updt_for <- function(pars, model, k, m) {
model$Q[1, 1, 1] <- exp(pars[1])
model$Q[2, 2, 1] <- exp(pars[2])
if (m == 1)
model$Q[3, 3, 1] <- exp(pars[3])
diag(model$Q[(3+m):(2+(k*2)+m),(3+m):(2+(k*2)+m), 1]) <- exp(pars[3+m])
model$H[1,1,1] <- exp(pars[4+m])
model
}
i <- 1
for (m in 0:1) {
for(k in harm){
if (m == 0) {
model <- SSModel(y_ucm ~ SSMtrend(1, NA) + SSMseasonal(7, NA, "dummy") + SSMseasonal(365, NA, "trig", harmonics = 1:k), H = NA)
name <- paste0("RW + seas_7 + seas_365 - k = ",k)
init <- numeric(4)
}
else {
model <- SSModel(y_ucm ~ SSMtrend(2, list(0,NA)) + SSMseasonal(7, NA, "dummy") + SSMseasonal(365, NA, "trig", harmonics = 1:k), H = NA)
init <- numeric(5)
init[1] <- 0
name <- paste0("IRW + seas_7 + seas_365 - k = ",k)
}
model$P1inf <- model$P1inf * 0 # trasformata in matrice di zero
model$a1[1] <- mean(y_ucm, na.rm = TRUE) # media del livello (media più sensata)
diag(model$P1) <- vary # a tutta la diagonale mettiamo la varianza della serie storica (upper bound)
init[(1+m)] <- log(vary) # log-var(dist.eta)
init[(2+m)] <- log(vary/10) # log-var(dist.seas) dummy
init[(3+m)] <- log(vary/10) # log-var(dist.seas) trig
init[(4+m)] <- log(vary) # log-var(err.oss)
fit <- fitSSM(model, init, updt_for, update_args = c(k = k, m = m))
smo <- KFS(fit$model, smoothing = c("state","mean"))
mape_train <- mean(abs(smo$muhat[1:(length(y_train))] - y_train_num)/y_train_num) * 100
mape_val <- mean(abs(smo$muhat[(length(y_train)+1):(nrow(dt))] - y_val_num)/y_val_num) * 100
best_ucm[i,] <- c(name, mape_train, mape_val)
i <- i + 1
}
}
write.csv(best_ucm,"csv_weights/best_ucm.csv", row.names = FALSE)
best_ucm <- read.csv("csv_weights/best_ucm.csv")
best_ucm[which.min(best_ucm$Mape_val),]
## Name Mape_train Mape_val
## 5 RW + seas_7 + seas_365 - k = 14 7.299592 15.18734
mod2_ucm <- SSModel(y_ucm ~ SSMtrend(1, NA) +
SSMseasonal(7, NA, "dummy") +
SSMseasonal(365, NA, "trig", harmonics = 1:14), H = NA)
init <- numeric(4)
init[1] <- log(vary) # log-var(dist.eta)
init[2] <- log(vary/10) # log-var(dist.seas) dummy
init[3] <- log(vary/10) # log-var(dist.seas) trig
init[4] <- log(vary) # log-var(err.oss)
mod2_ucm$P1inf <- mod2_ucm$P1inf * 0
mod2_ucm$a1[1] <- mean(y_ucm, na.rm = TRUE)
diag(mod2_ucm$P1) <- vary
fit2 <- fitSSM(mod2_ucm, init, updt_for, update_args = c(k = 14, m = 0))
#fit2$optim.out$convergence
smo2 <- KFS(fit2$model, smoothing = c("state","mean"))
print(paste0("MAPE on training set ", mean(abs(smo2$muhat[1:(length(y_train))] - y_train_num)/y_train_num) * 100))
## [1] "MAPE on training set 7.29959174329699"
print(paste0("MAPE on validation set ", mean(abs(smo2$muhat[(length(y_train)+1):(nrow(dt))] - y_val_num)/y_val_num) * 100))
## [1] "MAPE on validation set 15.1873440009321"
# plot validation
plot(y_val_num, type = "l")
lines(smo2$muhat, col = "blue")
Proviamo ad inserire nel modello le dummy per le festivitÃ
dummy_holidays_ucm <- dm[holidays][1:length(y_ucm),]
mod3_ucm <- SSModel(y_ucm ~ dummy_holidays_ucm$EasterSat +
dummy_holidays_ucm$EasterMon +
dummy_holidays_ucm$EasterTue +
dummy_holidays_ucm$Dec24 +
dummy_holidays_ucm$Dec25 +
dummy_holidays_ucm$Dec26 +
dummy_holidays_ucm$Dec31 +
dummy_holidays_ucm$Dec8 +
dummy_holidays_ucm$Jan1 +
dummy_holidays_ucm$Nov1 +
dummy_holidays_ucm$Jan6 +
dummy_holidays_ucm$Aug15 +
dummy_holidays_ucm$Apr25 +
dummy_holidays_ucm$Mag1 +
dummy_holidays_ucm$Giu2 +
SSMtrend(1, NA) + SSMseasonal(7, NA, "dummy") + SSMseasonal(365, NA, "trig", harmonics = 1:14), H = NA)
init <- numeric(4)
init[1] <- log(vary) # log-var(dist.eta)
init[2] <- log(vary/10) # log-var(dist.seas) dummy
init[3] <- log(vary/10) # log-var(dist.seas) trig
init[4] <- log(vary) # log-var(err.oss)
mod3_ucm$P1inf <- mod3_ucm$P1inf * 0
mod3_ucm$a1[1] <- mean(y_ucm, na.rm = TRUE)
diag(mod3_ucm$P1) <- vary
fit3 <- fitSSM(mod3_ucm, init, updt_for, update_args = c(k = 14, m = 0))
smo3 <- KFS(fit3$model, smoothing = c("state","mean"))
print(paste0("MAPE on training set ", mean(abs(smo3$muhat[1:(length(y_train))] - y_train_num)/y_train_num) * 100))
## [1] "MAPE on training set 6.51449913711938"
print(paste0("MAPE on validation set ", mean(abs(smo3$muhat[(length(y_train)+1):(nrow(dt))] - y_val_num)/y_val_num) * 100))
## [1] "MAPE on validation set 20.9862678179542"
# plot validation
plot(y_val_num, type = "l")
lines(smo3$muhat, col = "blue")
L’aggiunta dei regressori peggiora le performance del modello. Come si può vedere dal Mape del training set, il modello overfitta.
import numpy as np
import pandas as pd
import math
from keras.models import Sequential
## Using TensorFlow backend.
from keras.layers import Dense
from keras.layers import LSTM, GRU, Dropout
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error
def create_window(dataset, look_back):
X = list()
Y = list()
for i in range((len(dataset)-look_back)):
X.append(dataset[i:i+look_back])
Y.append(dataset[i+look_back])
return np.array(X), np.array(Y)
def GRU_model(look_back):
model = Sequential()
model.add(GRU(64, dropout = 0.1, recurrent_dropout = 0.5, input_shape=(1, look_back), return_sequences=True))
model.add(GRU(128, dropout = 0.1, recurrent_dropout = 0.5,activation="relu"))
model.add(Dense(1))
print(model.summary())
model.compile(loss='mean_squared_error', optimizer='Adam')
return model
dataset = r.y_box
dataset = dataset.astype('float32')
# normalize the dataset
scaler = MinMaxScaler(feature_range=(0, 1))
dataset = np.array(dataset)
dataset_scaled = scaler.fit_transform(dataset.reshape(-1,1))
train_size = len(dataset) - 365
val_size = len(dataset) - train_size
train_scaled, val_scaled = dataset_scaled[0:train_size,:], dataset_scaled[train_size:len(dataset_scaled),:]
batch_size = 1
look_back = 365
trainX_scaled, trainY_scaled = create_window(train_scaled, look_back)
trainX_scaled = np.reshape(trainX_scaled, (trainX_scaled.shape[0], 1, trainX_scaled.shape[1]))
model = GRU_model(look_back)
#model.fit(trainX_scaled, trainY_scaled, epochs=20, batch_size=1, verbose=1)
## Model: "sequential_1"
## _________________________________________________________________
## Layer (type) Output Shape Param #
## =================================================================
## gru_1 (GRU) (None, 1, 64) 82560
## _________________________________________________________________
## gru_2 (GRU) (None, 128) 74112
## _________________________________________________________________
## dense_1 (Dense) (None, 1) 129
## =================================================================
## Total params: 156,801
## Trainable params: 156,801
## Non-trainable params: 0
## _________________________________________________________________
## None
model.load_weights("csv_weights/model_GRU_stacked.h5")
def mean_absolute_percentage_error(y_true, y_pred):
y_true, y_pred = np.array(y_true), np.array(y_pred)
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
trainPredict_scaled = model.predict(trainX_scaled, batch_size=batch_size)
trainPredict = scaler.inverse_transform(trainPredict_scaled)
trainY = scaler.inverse_transform(trainY_scaled)
trainScore = mean_absolute_percentage_error(trainY, trainPredict)
print('Train Score: %.2f MAPE' % (trainScore))
## Train Score: 9.60 MAPE
valX = trainX_scaled[len(trainX_scaled)-1]
valX = np.append(valX[0], trainY_scaled[len(trainY)-1][0])
val_pred = []
for i in range(365):
valX_i = valX[-365:].reshape(1,1,-1)
prediction = model.predict(valX_i, batch_size=batch_size)
valX = np.append(valX_i, prediction)
val_pred.append(prediction)
valPredict = scaler.inverse_transform(np.array(val_pred).reshape(-1,1))
Yval = r.y_val
valScore = mean_absolute_percentage_error(Yval, valPredict)
print('Val Score: %.2f MAPE' % (valScore))
## Val Score: 10.87 MAPE
# plot validation
plot(y_val_num, type = "l")
lines(as.numeric(py$valPredict), col = "blue")
# plot training
plot(y_train_num[366:(length(y_train_num))], type = "l")
lines(as.numeric(py$trainPredict), col = "blue")
fs_train_final <- fourier(ts(y_box, frequency = 365), K = 1)
fs_test_final <- fourier(ts(y_box, frequency = 365), K = 1, h = 334)
dummy_holidays_train_final <- dm[holidays][1:length(y),]
dummy_holidays_test_final <- dm[holidays][(length(y)+1):(length(y)+334),]
mod5_reg_final <- Arima(y_box, c(6,1,1), list(order=c(1,1,1), period=7),
xreg=as.matrix(cbind(fs_train_final,dummy_holidays_train_final)), include.constant = TRUE)
#mod5_reg
#checkresiduals(mod5_reg)
pred5_reg_final <- forecast(mod5_reg_final, h=334, xreg=as.matrix(cbind(fs_test_final,dummy_holidays_test_final)))
print(paste0("MAPE on training set ", mean(abs(mod5_reg_final$fitted - as.numeric(y_box))/as.numeric(y_box)) * 100))
## [1] "MAPE on training set 8.37480241207164"
plot(as.numeric(pred5_reg_final$mean), type = "l")
#autoplot(pred5_reg_final)
y_ucm_final <- as.numeric(y_box)
y_ucm_final[(length(y_ucm_final)+1):(length(y_ucm_final)+334)] <- NA
updt_for <- function(pars, model, k, m) {
model$Q[1, 1, 1] <- exp(pars[1])
model$Q[2, 2, 1] <- exp(pars[2])
if (m == 1)
model$Q[3, 3, 1] <- exp(pars[3])
diag(model$Q[(3+m):(2+(k*2)+m),(3+m):(2+(k*2)+m), 1]) <- exp(pars[3+m])
model$H[1,1,1] <- exp(pars[4+m])
model
}
mod_ucm_final <- SSModel(y_ucm_final ~ SSMtrend(1, NA) + SSMseasonal(7, NA, "dummy") + SSMseasonal(365, NA, "trig", harmonics = 1:14), H = NA)
vary <- var(y_ucm_final, na.rm = TRUE)
init <- numeric(4)
init[1] <- log(vary) # log-var(dist.eta)
init[2] <- log(vary/10) # log-var(dist.seas) dummy
init[3] <- log(vary/10) # log-var(dist.seas) trig
init[4] <- log(vary) # log-var(err.oss)
mod_ucm_final$P1inf <- mod_ucm_final$P1inf * 0
mod_ucm_final$a1[1] <- mean(y_ucm_final, na.rm = TRUE)
diag(mod_ucm_final$P1) <- vary
fit_final <- fitSSM(mod_ucm_final, init, updt_for, update_args = c(k = 14, m = 0))
smo_final <- KFS(fit_final$model, smoothing = c("state","mean"))
print(paste0("MAPE on training set ", mean(abs(smo_final$muhat[1:(length(y))] - as.numeric(y_box))/as.numeric(y_box)) * 100))
## [1] "MAPE on training set 6.83152740363415"
plot(smo_final$muhat[(length(y_box)+1):(length(y_ucm_final))], type = "l",col = "blue")
validation_scaled = dataset_scaled[(train_size-365):len(dataset),:]
validationX_scaled, validationY_scaled = create_window(validation_scaled, look_back)
validationX_scaled = np.reshape(validationX_scaled, (validationX_scaled.shape[0], 1, validationX_scaled.shape[1]))
model_final = GRU_model(look_back)
## Model: "sequential_2"
## _________________________________________________________________
## Layer (type) Output Shape Param #
## =================================================================
## gru_3 (GRU) (None, 1, 64) 82560
## _________________________________________________________________
## gru_4 (GRU) (None, 128) 74112
## _________________________________________________________________
## dense_2 (Dense) (None, 1) 129
## =================================================================
## Total params: 156,801
## Trainable params: 156,801
## Non-trainable params: 0
## _________________________________________________________________
## None
model_final.load_weights("csv_weights/model_GRU_stacked_valid.h5")
#model_final.load_weights("model_GRU_stacked.h5")
#model_final.fit(validationX_scaled, validationY_scaled, epochs=4, batch_size=1, verbose=0)
# make predictions
trainPredict_final_scaled = model_final.predict(trainX_scaled, batch_size=batch_size)
trainPredict_final = scaler.inverse_transform(trainPredict_final_scaled)
# invert predictions
trainY = scaler.inverse_transform(trainY_scaled)
trainScore_final = mean_absolute_percentage_error(trainY, trainPredict_final)
print('Train Score: %.2f MAPE' % (trainScore_final))
## Train Score: 13.40 MAPE
# make predictions
validPredict_scaled = model_final.predict(validationX_scaled, batch_size=batch_size)
validPredict_final = scaler.inverse_transform(validPredict_scaled)
# invert predictions
# i dati di validation sono stati passati al modello per il training finale
validationY = scaler.inverse_transform(validationY_scaled)
validScore_final = mean_absolute_percentage_error(validationY, validPredict_final)
print('Train Score: %.2f MAPE' % (validScore_final))
## Train Score: 7.27 MAPE
testX = validationX_scaled[len(validationX_scaled)-1]
testX = np.append(testX[0], trainY_scaled[len(validationY)-1][0])
test_pred = []
for i in range(334):
test_i = testX[-365:].reshape(1,1,-1)
prediction = model_final.predict(test_i, batch_size=batch_size)
testX = np.append(test_i, prediction)
test_pred.append(prediction)
testPredict = scaler.inverse_transform(np.array(test_pred).reshape(-1,1))
plot(py$testPredict, type = "l", col="red")
pred_arima <- as.numeric(pred5_reg_final$mean)
pred_ucm <- smo_final$muhat[(length(y)+1):(length(y_ucm_final))]
pred_gru <- as.numeric(py$testPredict)
pred_arima_inv <- InvBoxCox(pred_arima, lambda, biasadj = TRUE, fvar = var(pred_arima))
pred_ucm_inv <- InvBoxCox(pred_ucm, lambda, biasadj = TRUE, fvar = var(pred_ucm))
pred_gru_inv <- InvBoxCox(pred_gru, lambda, biasadj = TRUE, fvar = var(pred_gru))
test_data <- dm$Date[(length(y_train_num)+1):(length(y_train_num)+334)]
pred_test <- cbind(test_data,pred_arima_inv, pred_ucm_inv,pred_gru_inv)
colnames(pred_test) <- c("Data","ARIMA","UCM","GRU")
write.csv(pred_test, "csv_weights/SDMTSA_847160.csv",)